/*
  3D - C++ Class Library for 3-D Rigid Transformations
  Copyright (C) 1996, 1997  Gino van den Bergen

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Library General Public
  License as published by the Free Software Foundation; either
  version 2 of the License, or (at your option) any later version.

  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Library General Public License for more details.

  You should have received a copy of the GNU Library General Public
  License along with this library; if not, write to the Free
  Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

  Please send remarks, questions and bug reports to gino@win.tue.nl,
  or write to:
                  Gino van den Bergen
		  Department of Mathematics and Computing Science
		  Eindhoven University of Technology
		  P.O. Box 513, 5600 MB Eindhoven, The Netherlands
*/

#ifndef _VECTOR_H_
#define _VECTOR_H_

#include <assert.h>
#include <iostream.h>

#include "Basic.h"

class Vector {
public:
  Vector() {}
  Vector(const float v[3]) { setValue(v); }
  Vector(const double v[3]) { setValue(v); }
  Vector(Scalar x, Scalar y, Scalar z) { setValue(x, y, z); }
  
  Scalar  operator[](int i) const { return coor[i]; }
  Scalar& operator[](int i)       { return coor[i]; }
  
  void setValue(const float v[3]) {
    coor[X] = v[X]; coor[Y] = v[Y]; coor[Z] = v[Z];
  }
  
  void setValue(const double v[3]) {
    coor[X] = v[X]; coor[Y] = v[Y]; coor[Z] = v[Z];
  }
  
  void setValue(Scalar x, Scalar y, Scalar z) {
    coor[X] = x; coor[Y] = y; coor[Z] = z;
  }
  
  operator const Scalar *() const { return coor; }
  
  Vector& operator+=(const Vector& v);
  Vector& operator-=(const Vector& v);
  Vector& operator*=(Scalar s);
  Vector& operator/=(Scalar s);
  
  bool approxZero() const;

  int closestAxis() const;

  Scalar length2() const;
  Scalar length() const;

  void normalize();
  Vector normalized() const;

  static Vector random();

protected:
  Scalar coor[3];                            
};

Vector operator+(const Vector& v1, const Vector& v2);
Vector operator-(const Vector& v1, const Vector& v2);
Vector operator-(const Vector& v);
Vector operator*(const Vector& v, Scalar s);
Vector operator*(Scalar s, const Vector& v);
Vector operator/(const Vector& v, Scalar s);

Scalar dot(const Vector& v1, const Vector& v2);

Scalar length2(const Vector& v);
Scalar length(const Vector& v);

Scalar angle(const Vector& v1, const Vector& v2);
Vector cross(const Vector& v1, const Vector& v2);

ostream& operator<<(ostream& os, const Vector& v);




inline Vector& Vector::operator+=(const Vector& v) {
  coor[X] += v[X]; coor[Y] += v[Y]; coor[Z] += v[Z];
  return *this;
}

inline Vector& Vector::operator-=(const Vector& v) {
  coor[X] -= v[X]; coor[Y] -= v[Y]; coor[Z] -= v[Z];
  return *this;
}
 
inline Vector& Vector::operator*=(Scalar s) {
  coor[X] *= s; coor[Y] *= s; coor[Z] *= s;
  return *this;
}

inline Vector& Vector::operator/=(Scalar s) {
  assert(!eqz(s));
  coor[X] /= s; coor[Y] /= s; coor[Z] /= s;
  return *this;
}

inline Vector operator+(const Vector& v1, const Vector& v2) {
  return Vector(v1[X] + v2[X], v1[Y] + v2[Y], v1[Z] + v2[Z]);
}

inline Vector operator-(const Vector& v1, const Vector& v2) {
  return Vector(v1[X] - v2[X], v1[Y] - v2[Y], v1[Z] - v2[Z]);
}

inline Vector operator-(const Vector& v) {
  return Vector(-v[X], -v[Y], -v[Z]);
}

inline Vector operator*(const Vector& v, Scalar s) {
  return Vector(v[X] * s, v[Y] * s, v[Z] * s);
}

inline Vector operator*(Scalar s, const Vector& v) { return v * s; }

inline Vector operator/(const Vector& v, Scalar s) {
  assert(!eqz(s));
  return Vector(v[X] / s, v[Y] / s, v[Z] / s);
}

inline bool Vector::approxZero() const {
  return eqz(coor[X]) && eqz(coor[Y]) && eqz(coor[Z]);
}

inline int Vector::closestAxis() const {
  Scalar a[2];
  a[X] = abs(coor[X]);
  a[Y] = abs(coor[Y]);
  int axis = a[X] < a[Y] ? Y : X;
  return a[axis] < abs(coor[Z]) ? Z : axis;
}

inline Scalar dot(const Vector& v1, const Vector& v2) {
  return v1[X] * v2[X] + v1[Y] * v2[Y] + v1[Z] * v2[Z];
}

inline Scalar Vector::length2() const { return dot(*this, *this); }
inline Scalar Vector::length() const { return sqrt(length2()); }

inline void Vector::normalize() { *this /= length(); }
inline Vector Vector::normalized() const { return *this / length(); }

inline Vector Vector::random() {
  Scalar z = 2 * rnd() - 1;
  Scalar r = sqrt(1 - z * z);
  Scalar t = TWO_TIMES_PI * rnd();
  return Vector(r * cos(t), r * sin(t), z);
}

inline Scalar length2(const Vector& v) { return v.length2(); }
inline Scalar length(const Vector& v) { return v.length(); }

inline Scalar angle(const Vector& v1, const Vector& v2) {
  Scalar s = sqrt(v1.length2() * v2.length2());
  assert(!eqz(s));
  return acos(dot(v1, v2) / s);
}

inline Vector cross(const Vector& v1, const Vector& v2) {
  return Vector(v1[Y] * v2[Z] - v1[Z] * v2[Y],
		v1[Z] * v2[X] - v1[X] * v2[Z],
		v1[X] * v2[Y] - v1[Y] * v2[X]);
}

inline ostream& operator<<(ostream& os, const Vector& v) {
  return os << v[X] << ' ' << v[Y] << ' ' << v[Z];
}

#endif
